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Abstract 

The aims of the study were to establish the prevalence of high bone mass (HBM) in a cohort of Spanish postmenopausal 
women (BARCOS) and to assess the contribution of LRP5 and DKKl mutations and of common bone mineral density (BMD) 
variants to a HBIVl phenotype. Furthermore, we describe the expression of several osteoblast-specific and Wnt-pathway 
genes in primary osteoblasts from two HBM cases. A 0.6% of individuals (10/1600) displayed Z-scores in the HBM range (sum 
Z-score >4). While no mutation in the relevant exons of LRP5 was detected, a rare missense change in DKKl was found 
(p.Y74F), which cosegregated with the phenotype in a small pedigree. Fifty-five BMD SNPs from Estrada ef al. [NatGenet 
44:491-501,2012] were genotyped in the HBM cases to obtain risk scores for each individual. In this small group of samples, 
Z-scores were found inversely related to risk scores, suggestive of a polygenic etiology. There was a single exception, which 
may be explained by a rare penetrant genetic variant, counterbalancing the additive effect of the risk alleles. The expression 
analysis in primary osteoblasts from two HBM cases and five controls suggested that IL6R, DLX3, TWISTl and PPARG are 
negatively related to Z-score. One HBM case presented with high levels of RUNX2, while the other displayed very low S0X6. 
In conclusion, we provide evidence of lack of LRP5 mutations and of a putative HBM-causing mutation in DKKl. Additionally, 
we present SNP genotyping and expression results that suggest additive effects of several genes for HBM. 
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Introduction 

Osteoporosis has a complex genetic background. Bone mineral 
density (BMD) is a highly heritable intermediate phenotype that 
correlates well with fracture risk [1-6]. BMD is distributed as a 
Gaussian curve in the general population, with two small groups 
having extremely low or extremely high BMD values at both ends. 
These individuals with extreme phenotypes may bear infrequent 
and highly penetrant alleles at a few specific loci. Alternatively, the 
extreme phenotypes may depend on the presence of a high 
number of common variants with low penetrance and additive 
effects. 

A few individuals with high bone mass (HBM, MIM#601884), 
as defined by a sum Z-score >4 (total lumbar spine Z-score -I- total 
femoral neck Z-score), have been reported to bear highly 
penetrant missense alleles at the low-density lipoprotein recep- 
tor-related protein 5 {LRP5, MIM#603506) locus that are 



transmitted in an autosomal dominant way. More than 10 years 
ago, two different groups found that LRP5 regulated bone mass 
[7,8]. While inactivating mutations in LRP5 were shown to cause 
osteoporosis-pseudoghoma syndrome [7], gain-of-function muta- 
tions caused a high bone mass (HBM) phenotype [8]. This 
phenotype has been associated with the LRP5-G 171V mutation in 
two independent pedigrees [8,9]. Six additional missense muta- 
tions (Dl 1 lY, G171R, A214T, A214V, A242T and T253I), all in 
the first P-propeUer domain of LRP5, were identified in patients 
who also showed an increased bone density [10]. The affected 
individuals had elevated bone synthesis assessed by serum markers, 
but normal bone resorption, bone architecture and serum calcium, 
phosphate, PTH and vitamin D levels [8,9]. Significant pheno- 
typic heterogeneity was reported, and some affected family 
members also had a torus palatinus. 
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LRP5 acts as a co-receptor with members of the Frizzled family 
to activate the canonical W nt/ P-catenin signalling pathway, which 
is crucial for bone formation [1 1]. This pathway is activated by the 
binding of the appropriate Wnt protein to LRP5 and is blocked by 
the binding of inhibitors such as Dickkopf-related protein 1 
(encoded by DKKl, MIM#605189) and Sclerostin (encoded by 
SOST). The HBM-causing mutation prevents the binding of these 
two inhibitors. Mutations in SOST are the cause of van Buchem 
disease [12] and sclerostosis [13], two pathologies with an 
abnormally high bone density. On the other hand, Dkkl'^'~ mice 
showed a marked increase in bone mass [14]. 

The prevalence of HBM in the general population has been 
estimated as 0.2-1% [8,15,16], but the genetic architecture of this 
extreme phenotype remains poorly understood. However, recent 
genome-wide association (GWA) analyses and meta-analyses have 
established a number of genomic loci that explain differences in 
BMD across the general population. In particular, Estrada et al. 
[1 7] identified 56 such genomic loci and showed how they can be 
used to calculate risk scores to predict BMD. 

In order to explore the genetic constitution of a high bone mass 
phenotype, our aims were, first, to establish the prevalence of 
HBM in the BARCOS (BARCelona Osteoporosis) cohort of 
postmenopausal Spanish women; second, to determine whether 
any of the HMB cases carried LRP5 or DKKl mutations that could 
explain the phenotype; and third, to assess whether the HBM cases 
were carriers of a low number of risk alleles of 55 autosomal 
GWA-identified BMD loci. Also, we took advantage of the 
availability of primary osteoblasts from two HBM cases to 
characterize the osteoblast RNA in terms of osteoblast-specific 
and/or Wnt-pathway genes by comparison with osteoblast RNA 
from donors with normal or low BMD. 

Materials and Methods 

Ethics Statement 

Both the Bioethics Committee of Universitat de Barcelona and 
the Clinical Research Ethics Committee of Pare de Salud MAR 
have emitted a favourable bioethical statement regarding the 
present research. Specifically, the protocol for obtention of 
peripheral blood from the BARCOS cohort women and the 
protocol for the obtention of primary osteoblasts from bone 
specimens extracted from knee samples otherwise discarded at the 
time of orthopaedic surgery, were approved by both committees. 
Written informed consents were obtained from the participants in 
both instances. 

Study Cohort 

The study population (listed in Table 1) included the HBM 
cases in the BARCOS cohort (n = 10 unrelated cases). This cohort 
of postmenopausal women from the Barcelona area has been 
described elsewhere [18,19]. At present time, it includes DXA 
values for 1 600 women and DNA samples for 1001 of them. Six 
additional unrelated female HBM cases were recruited from 3 
hospitals of Barcelona, Hospital de Sant Pau (1 case), CETIR (a 
private medical services centre specializing in nuclear medicine 
and other imaging modalities, 4 cases), and Hospital de 
L'Esperan^a (1 case). Some relatives of particular cases were also 
studied. Blood samples and written informed consent were 
obtained in accordance with the regulations of the Clinical 
Research Ethics Committee of Pare de Salud MAR. A total of 
1600 dual-energy X-ray absorptiometry scans (DXA; QDR 4500 
SL; Hologic, Waltham, MA, USA) of the women from this cohort 
were analysed in order to pinpoint those HBM cases in which the 
sum Z-score (hip plus lumbar spine) was equal to or greater than 



four [8]. All DXA measurements were performed prior to any 
treatment that could increase bone mass. Pathologic phenotypes 
such as osteopetrosis or any other sclerosing bone disorders were 
ruled out based on the absence of radiologic alterations in the skuU 
and long bones, the absence of any fragility fracture and the 
absence of any underlying disease. For individuals HBMl, HBM8, 
and HBMl 1, no DNA sample was available. 

Sequencing of LRP5 and DKKl 

The genomic DNA of n = 13 HBM cases (Table 2) was isolated 
from peripheral blood leukocytes using conventional methods. In 
all probands, LRP5 exons 2 to 4 (encoding the first beta-propeUer 
and harbouring all HBM-related mutations described so far) and 
exons 9 to 16 (encoding the third and fourth beta-propeUers 
involved in binding to DKKl), and their intronic flanking regions, 
and the four exons and flanking regions of DKKl were amplified 
and sequenced using specific primers. Mutation screening was 
performed by direct sequencing using the BigDye v3. 1 kit (Applied 
Biosystems, Foster City, CA, USA) and the ABI PRISM 3730 
DNA Analyzer (Applied Biosystems). Primers and PCR conditions 
are listed in Table SI. Nomenclature for DNA variants followed 
these reference sequences: NM_002335.2 for LRP5 and 
NM_012242.2 ior DKKl. 

SNP Selection and Genotyping 

Out of the 64 SNPs identified by Estrada et al. [17], we chose 
55, one for each autosomal locus (listed in Table S2). Genotyping 
of the 13 available HBM cases was carried out with a KASPar 
v4.0 genotyping system at the Kbioscience facilities (KBioscience, 
Herts, UK) using the Kraken aUele-caUing algorithm [20]. The 
genotypes of 1001 BARCOS participants for the same SNPs were 
already available, since this cohort was included in the replication 
phase of the study by Estrada et al. [17]. One of the SNPs 
(rs3790160) gave conflicting results and was eliminated from the 
analyses. Quality control was carried out by resequencing 6.28% 



Table 1 

women. 


. Z-score values, age 


and CO 


hort of the participant 




Case 


Surr 


Z-score LS Z-score HIP Z-score Age 


Cohort 


hbmV- ^ 


6.1 


3.4 


2.7 


51 


BARCOS 


HBM2 


4.6 


3.0 


1.6 


57 


BARCOS 


HBM3 


4.9 


2.5 


2.4 


55 


BARCOS 


HBM4 


4.5 


2.5 


2.0 


62 


BARCOS 


HBM5 


4.5 


2.4 


2.1 


66 


BARCOS 


HBM6 


5.1 


2.6 


2.5 


52 


BARCOS 


HBM7 


4.6 


2.5 


2.1 


61 


BARCOS 


HBMS^' 


7.9 


4.0 


3.9 


55 


BARCOS 


HBM9 


7.0 


4.6 


2.4 


66 


BARCOS 


HBM10 


5.1 


2.8 


2.3 


75 


BARCOS 


HBMll^ 


6.8 


3.8 


3.0 


55 


HSANTPAU 


HBM12 


6.4 


3.8 


2.6 


59 


CETIR 


HBM13 


5.2 


2.6 


2.6 


67 


CETIR 


HBM14 


6.0 


3.7 


2.3 


64 


CETIR 


HBM15 


4.5 


2.4 


2.1 


54 


CETIR 


HBM16 


5.3 


3.6 


1.7 


77 


HESP 



^Deceased during the course of the study. 

^No DNA sample available. 

doi:l 0.1 371 /journal.pone.0094607.t001 



PLOS ONE I www.plosone.org 



2 



April 2014 I Volume 9 | Issue 4 | e94607 



Genetic Analysis of High Bone IVIass 



Table 2. Genotypes at LRP5 and DKKl. 









LRP5 exonic SNPs 






DKKl exonic 
changes 


Case 


P.V667M' 


P.A1330V' 


P.V1119V' 


Others' 


P.A106A' Variant 


HBM2 






Het. 


P.E644E 


Het. 


HBM3 






Het. 


P.N705N 


Homo. 


HBM4 




Het. 


Het. 


P.N740N 


Het. 


HBM5 






Het. 


P.E644E 




HBM6 










Homo. 


HBM7 






Het. 


P.N740N 




HBM9 




Het. 


Het. 




Het. 


HBMIO 






Het. 






HBM12 






Het. 


P.E644E 


Het. 


HBM13 


Het. 


Het. 


Het. 


P.N740N 


Het. 


HBM14 




Het. 


Het. 


P.N740N 


Het. 


HBM15 










p.Y74F^ 


HBM16 








P.E644E 







^Corresponding reference sequences, rs-numbers and MAFs are: LRP5 
(NM_002335.2): p.E644E {rs2277268. 0.06); p.V667 (rs4988321. 0.03); p.N705N 
(rsl 45456776. <0.01); p.N740N {rs2306862. 0.15); p.VlllW {rs556442. 0.28); 
P.A1330V (rs3736228. 0.13). DKKl (NM_01 2242.2): p.A106A (rs2241529. 0.46). All 
LRP5 variants listed under "Others", as well as the DKKl p.Y74F were found in 
heterozygous state. 

^A novel missense change is indicated in bold. Het.: heterozygous for the 
variant; Homo.: homozygous for the minor allele; NA: not available; -: 
homozygous for the reference allele. 
doi:l 0.1 371 /journal.pone.0094607.t002 

of the samples. The readings showed full concordance between the 
two techniques. 

Genetic Risk Allele Analysis 

Fifty-five SNPs previously described to be associated with BMD 
at GWA significance [17] were geno typed in the 13 HBM cases 
and the genetic risk score for each individual was then calculated 
by taking into account both the number of risk alleles and the 
relative effect of each SNP on BMD, as carried out by Estrada et 
al. [17]. This calculation was also performed for 1001 individuals 
in the BARCOS cohort with available SNP genotype and LS- 
BMD data. Briefly, the genotype of each SNP was transformed 
into a risk score by taking into account the effects estimated by the 
authors and listed in their Supplementary Table 9. The effect size 
(beta parameter or slope) is the BMD decrease due to the presence 
of one copy of the risk allele. The scores for homozygotes for the 
risk allele were 2x the effect size; scores for heterozygotes were Ix 
the effect size; and the scores for homozygotes for the alternative 
allele were zero. For each individual, the risk scores of all SNPs 
were summed up to obtain a global risk score, which was then 
normalized by dividing it by the mean effect size of BARCOS. 
Normalized global risk scores were sorted into five bins, as 
described in Estrada et al. [17]. Missing genotypes within 
BARCOS cohort individuals and also within HBM probands 
were solved by replacing them by the mean of the corresponding 
SNP scores in BARCOS. This strategy would attenuate the 
variance of the overall group [21]. 

Primary Osteoblast Isolation and Cell Culture 

Primary osteoblast (hOB) cells of postmenopausal women were 
available from bone specimens extracted from knee samples that 



would otherwise have been discarded at the time of artroplasty. 
Both informed consent and BMD values were obtained from the 
donors. Two of them were patients HBMIO and HBM16. Five 
female donors with sum Z-score values (and ages) of 2.4 (85 y), 1.2 
(74 y), 0.4 (79 y), -0.7 (70 y) and -2.2 (61 y), were used as 
controls. The bone tissue used was obtained from a region at least 
2 cm apart from the subchondral and the osteochondral plates far 
from the described 6 mm layer of trabecular bone below the above 
mentioned plates [22]. With this we aimed to minimize the 
potential issues of using osteoblastic cells located close to the 
damaged joint (given the known relationship and mutual 
influences between the inflamed cartilage and adjacent bone cells) 
and the alterations in phenotype and gene expression of the 
osteoblasts in this zone [23,24]. The hOB cells were obtained, as 
described previously [25,26]. Briefly, the trabecular bone was 
separated and cut into small fragments, washed in phosphate 
buffered solution (PBS) to remove non-adherent cells, and placed 
on a petri dish. Samples were incubated in Dulbecco's Modified 
Eagle Medium (DMEM; Gibco; Invitrogen, Paisley, UK), 
supplemented with sodium pyruvate (1 mM), L-glutamine 
(1 mM), 1% penicillin/streptomycin, 10% fetal calf serum (FCS), 
0.4% fungizone and 1% ascorbic acid. This allowed osteoblastic 
precursor cells to migrate from the fragments and proliferate. After 
confluence, cells were trypsinized and cultured in the same 
medium. When sub-confluence was reached again, the medium 
was aspirated and fresh medium with 10% serum or 0.1% bovine 
serum albumin was added. Forty-eight hours later, the medium 
was aspirated, cells were rinsed with phosphate-buffered saline, 
and the RNA was extracted. Standard histochemistry or 
quantitative PCR tests to measure alkaline phosphatase or 
osteocalcin expression, respectively, were used to confirm the 
osteoblastic nature of the cells, as described in [27]. 

RNA Extraction, cDNA Synthesis and Real-time PCR 

Total cell RNA was extracted using the High Pure RNA 
Isolation Kit (Roche Diagnostics, Mannheim, Germany) in 
accordance with the manufacturer's instructions. Two micrograms 
of total RNA were reverse transcribed using random primers of 
the High Capacity cDNA RT Kit (Applied Biosystems) in 
accordance with the manufacturer's instructions. The real-time 
PCR (qPCR) reactions were performed in a final volume of 10 Jil 
using 20 ng of each cDNA, which was used as template for each 
well in the RealTime ready Custom Panel 384 (Roche Diagnos- 
tics). This custom panel included 88 genes selected by us, taking 
into account, among other sources, the sites recently highlighted 
by GWA analyses, in particular those in Estrada et al. [17] and 
Duncan et al. [28]. All qPCR reactions for each sample were 
performed in triplicate with the LightCycler 480 Real-Time PCR 
System (Roche Diagnostics). Beta-2-microglobulin [B2M) was 
chosen as the reference gene because of its minimum coefficient 
of variation between samples. 

Validation of the 11 genes with positive results in the expression 
analysis described previously was performed with new assays 
designed using the online ProbeFinder software (Roche Diagnos- 
tics). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and 
18 S ribosomal {18 S) were tested as possible reference genes, and 
GAPDH was selected. For this validation step, samples from three 
additional control individuals were included in the analysis. Again, 
all qPCR reactions for each sample were performed in triplicate 
with the LightCycler 480 Real-Time PCR System (Roche 
Diagnostics). 
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Statistical Analysis 

Linear regression was analysed by first testing the most 
important assumptions (normality and homoscedasticity). Calcu- 
lations were performed with SPSS vll.5 (SPSS, Chicago, IL, 
USA). To test skewness (or asymmetry of a distribution) of the 
genetic risk distributions of the HBM cases and the BARCOS 
controls, we used the robust medcouple (MC) measure, with left 
and right tail weight measures (LMC and RMC), and constructed 
the MC-LR confidence interval [29]. 

Results 

HBM Prevalence in the BARCOS Cohort and Features of 
the HBIVl Cases 

In total, 1600 DXA scans were analysed across the BARCOS 
cohort. Those cases in which the sum Z-score was equal to or 
greater than four were considered HBM cases and further 
analysed. Pathologic phenotypes were ruled out based on a more 
in-depth examination of the medical history, a physical examina- 
tion and a radiologic study. In the BARCOS cohort, 10 cases 
(0.63% of individuals) fulfilled this HBM criterion. Six additional 
HBM cases were recruited elsewhere (see Materials and Methods). 
Z-score values, age and cohort for all HBM cases are listed in 
Table 1. 

Search for Mutations in LRP5 and DKK1 

Exons 2 to 4 of LRP5, which encode the first beta-propeller of 
the protein, and in which HBM mutations have pre\ iously been 
described, were sequenced in the 13 HBM cases with available 
DNA sample. Next, exons 9 to 16 were analysed, since they 
encode beta-propeUers 3 and 4, which have been described as 
binding regions for the LRP5 inhibitor DKKl. No novel or 
previously described causing mutations were found in any of these 
exons. The missense ZilPJ polymorphisms p.V667M (rs4988321) 
and P.A1330V (rs3736228), associated with BMD in GWA 
studies, together with other silent exonic variants found in the 
HBM individuals, are shown in Table 2. Their frequencies in 
HBM cases were similar to those found in the general population 
(dbSNP; http://www.ncbi.nlm.nih.gov/SNP/). Regarding DKKl, 
the four exons and flanking regions were amplified and sequenced 
in the 13 HBM cases. One previously undescribed heterozygous 
missense c:hange (p. Y74F) and an exonic silent polymorphism were 
found in different indi\iduals (Table 2). The p.Y74F was not 
present in the 1000 genomes database, while the tyrosine 74 and 
adjacent residues were found to be conserved in the Dkkl 
sequences of primates, rodents and cows, but not in C. lupus or D. 
rem (Fig. lA). PolyPhen-2 (HumDiv) and SIFT scores for this 
missense change were as follows: 0.48 (possibly damaging) and 
0.38 (tolerated), respectively. Examination of the offspring of case 
HBM 1.5 revealed a cosegregation of this mutation with the HBM 
phenotype (Fig. IB). The daughter was found to be an HBM case 
(sum Z-score: 4.9) and was heterozygous for the mutation, while 
her brother had a normal sum Z-score value (0.5) and did not 
carry this DNA change. 

Analysis of 55 Bone Mineral Density Loci 

Figure 2A shows the distribution of the BARCOS individuals 
into five different osteoporosis risk score bins (bars) and the mean 
LS-BMD value for each bin (triangles). Risk scores were derived 
from the number and effect of the BMD-associated SNPs 
described in Estrada et al. [17] (see Methods). As expected, a 
decrease in BMD values was observed as the genetic risk score 
increased (Pearson correlation coefficient =—0.972, p = 0.0057, 
r^ = 0.94). Figure 2B shows a similar graph for the 1 1 HBM 



individuals for which genotyping was successful (two of the HBM 
cases -HBM 14 and HBM 16- had to be discarded because of sub- 
optimal genotyping results). The frequency distribution shows a 
shift towards the lower genetic risk score bins. Again, BMD values 
(measured as Z-scores) decreased as genetic risk scores increased, 
with the exception of one individual (HBM9), who presented the 
maximum BMD value (Z-score = 7) and the largest genetic risk 
score. Interestingly, the HBM phenotype in the family of HBM9 
seems to segregate as a discrete trait: the mother is also a HBM 
individual (sum Z-score = 4.4), while her eldest brother is not 
(Fig. 2C). A correlation analysis between risk scores and Z-scores 
of the 1 1 HBM individuals was not possible due to the small 
sample size, since linear regression assumptions were not fulfilled. 
To compare the two distributions their skewness was analyzed. 
MC-LR 95% confidence intervals for skewness [(-0.193, 0.116) 
for Fig 2A and (-1.528, 0.740) for Fig 2B] pointed to a loss of 
symmetry for Fig 2B. However the intervals did overlap, due to the 
limited sample size of the HBM group. Thus, significant 
difierences between the two groups could not be formally 
demonstrated. 

Expression Analysis of 88 Bone-development and/or 
Wnt-pathway Genes 

A transcriptomic analysis by qPCR was carried out in primary 
osteoblast samples from two HBM and two age- and gender- 
matched control individuals that were obtained after knee- 
replacement orthopaedic surgery. Because of the small sample 
size, the approach was only meant to be descriptive. In a first step, 
1 1 out of 88 bone-development and/ or Wnt-pathway genes were 
selected (Table 3) due to differences above or below 2-fold in mean 
expression level between the two HBM cases and the two control 
individuals. Subsequently, these 1 1 genes were re-analysed using 
samples from the two HBM and five control individuals (the initial 
2 controls, whose Z scores were negative, and 3 new ones with Z- 
scores of 0.4, 1.2 and 2.4). Because the sum Z-scores of the control 
individuals were scattered across a wide range of values, expression 
levels were plotted against sum Z-scores and we observed that 
TWIST, ILR6, DLX3 and PPARG displayed a trend of correlation 
(Figure 3,A-D). R^ values for the four genes were 0.676, 0.895, 
0.807 and 0.461, respectively. However, due to the small sample 
size, linear regression assumptions were not fulfilled. Thus, no p- 
values are provided. We also observed that for S0X6 and RUNX2, 
one of the HBM samples (but not the other) presented an 
expression level that was 5-fold decreased or increased, respec- 
tively, compared with controls (Fig. 3E and 3F, respectively). 

Discussion 

In this study we established that the prevalence of HBM in a 
Spanish cohort of postmenopausal women is 0.63%. None of the 
HBM indi\'iduals had mutations in the relevant exons of the LRP5 
gene that could explain their phenotype. One individual had a 
rare missense change in DKKl (p.Y74F). The results of the analysis 
of 55 osteoporosis-predisposing SNPs pointed to an inverse 
correlation between risk alleles and BMD in this group of HBM 
women, with the exception of one case with the highest BMD 
value and the highest risk score. Finally, the results of an 
expression analysis in primary osteoblasts showed a negative trend 
between Z-scores and mRNA levels of TWISTl, IL6R, DLX3wA 
PPARG. 

There are few studies that describe the prevalence of HBM in 
the general population. In a recent one [15], the authors studied a 
UK DXA-scanned population in which the prevalence of HBM 
was 0.2% of individuals. The lower prevalence (0.2% versus our 
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H 


sapiens 


60 


P 


troglodytes 


60 


M 


mulatta 


60 


B 


taurus 


60 


M 


musculus 


61 


R 


norvegicus 


61 


C 


lupus 


42 


D 


rerio 


44 



VSAAPGILYPGGNFYQTIDNYQPYPCAEDE 89 

VSAAPGILYPGGNEYQTIDNYQPYPCAEDE 89 

VSAAPGI LYPGGNE Y QT IDNYQPYPCAEDE 8 9 

VSAAPGIPFEGGSFYQTIDNHQPYPCAEDE 89 

VSVAPGVLYEGGNK Y QTLDNYQPYPCAEDE 9 0 

VSVAPGVLYEGGNEYQTLDNYQPYPCAEDE 90 

SYVRGTVTFGRGSWGLTARSPQPYPCAEDE 71 

VSPSPDVSPLDSLNFALDTPQQPLICESDE 73 



B 




P.Y74F/+ 



0.5 



P.Y74F/+ +/+ 



Figure 1. Mutation p.Y74F of DKK1 may be responsible for high bone mass in family HBM15. (A) Alignment of a partial human DKKl 
sequence with those of several vertebrates. The tyrosine-74 residue is boxed. (B) Pedigree of case HBIVI1 5 (arrow): filled symbols indicate a high bone 
mass phenotype. Numbers inside symbols are sum Z-score values. 
doi:1 0.1 371/journal.pone.0094607.g001 



0.6%) may be due to difFerences in the study design, including the 
definition of HBM, which was stricter in their study. All HBM- 
related LRPS mutations identified to date are located in the first 
beta-propeUer of the protein [10], while the mutations that cause 
ostcoporosis-pseudogiioma syndrome and exudati\'e \ itrcorrtino- 
pathy arc found all over the gene [30] . In our mutational analysis 
of exons coding for the first beta-propeller (and for the third and 
fourth, which have been described as binding regions for DKKl 
[31]) no mutations were detected in the 13 HBM individuals 
analysed. Our results are in agreement with those published by 
Duncan et al. [32], who pointed out that <2% of their HBM cases 
were due to mutations on exons 2-4 and intron/ exon boundaries 
of the LRP5 gene, after analysing 98 patients. 

We also analysed the DKKl gene under the hypothesis that loss- 
of-function mutations in this gene could have the same effect as 
gain-of-function mutations in LRPS. In this regard, it has been 
shown that bone mass was in\-cr,scly proportional to Dkkl lcv('ls in 
mice [33], and there are therapies under development based on 
DKKl inhibition to increase bone mass (reviewed in Ke et al. [34]). 
We found a missense change (p.Y74F) in heterozygosis in one 
HBM individual. We have gathered some evidence that supports a 
causative role for this mutation. The change affects a conserved 
residue and it is predicted to be possibly damaging. Additionally, it 
cosegregates with HBM in the nuclear pedigree of case HBM 15. 
In humans, the only report on DKKl mutations is by Korvala et al. 
[35], who recentiy suggested that a mutation in DKKl may 
predispose individuals to primary osteoporosis. No mutations were 



found in the HGMD Professional 2012.4 database (released 29 
March 2013), and a limited number of very rare missense changes 
(31) were found in the 1000 genomes database (released 13 
December 2012); eleven of these are predicted to be deleterious by 
SIFT and probably damaging by PolyPhen. Whether these 
changes, or the p.Y74F described in this paper, are associated 
with a HBM phenotype remains an interesting open question, 
which win require functional analyses to be confirmed. 

To test whether HBM could be explained by a polygenic 
additive effect of susceptibility loci, we chose to genotype the 
GWA-discovered BMD loci defined by Estrada et al. [17]. We 
were able to use the BARCOS cohort information for these same 
loci as a framework for comparison. We would expect HBM 
individuals to bear a high number of protective alleles or a low 
number of risk alleles. As seen in Figure 2, the distribution of 
genetic risk scores for the HBM cases shifted towards lower risk, 
with the mode bin set at 44-48 instead of 48-52 and lost 
symmetry. However, due to the limited number of HBM cases, 
statistical significance was not achieved and, thus, no significant 
differences between the small cohort of individuals with HBM and 
the BARCOS cohort could be demonstrated. 

Additionally, Z-score values decreased as genetic risk scores 
increased, suggesting that common variation is playing a role in 
determining HBM. However, we note that the only HBM 
individual falling into the highest risk score bin (HBM9) was the 
one with the highest Z-score and such a contradictory fact might 
indicate the existence of a highly penetrant and probably rare 
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Figure 2. Distribution of genetic risl< scores in BARCOS and in tlie IHBM group. Distributions of genetic risl< scores among 1001 BARCOS 
individuals (A) and 1 1 HBIVI probands (B), and their relationships with BMD or Z-score values, are shown. Histograms describe counts of individuals in 
each genetic score category (left axis scale); (A) From left to right, exact numbers of individuals in each bin are: 88, 222, 355, 238 and 98. Triangles 
(right axis scale) represent LS-BIVID means and vertical bars depict their standard errors; (B) Diamonds represent mean Z-score values. (C) Pedigree of 
family HBM9. Arrow indicates the proband HBM9; filled symbols represent presence of the HBM phenotype; numbers below symbols denote sum Z- 
scores; NA: not available. 
doi:1 0.1 371/journal.pone.0094607.g002 
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Table 3. Eleven genes selected from the RealTime Custom 
Paner. 





Wnt-pathway genes 


Bone biology genes 




FZD3 


BMP4 IL6R 


S0X6 


SOST 


C0LWA1 PPARG 


SP7 (OSX) 




DLX3 RUNX2 


TWIST! 



^Those displaying at least a two-fold difference between the mean expression levels 
of HBM10 and HBM16 and the mean of the two controls with negative Z-scores. 
doi:l 0.1 371 /journal.pone.0094607.t003 
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protective allele counterbalancing the additive effect of the risk 
alleles. The HBM phenotype of this individual might be a 
Mendelian trait due to an as-yet unidentified gene. Altogether, 
these results are consistent with the coexistence of both polygenic 
and Mendelian cases of HBM, which would then be a 
heterogeneous trait. 

To gain further insight into specific genes that might be involved 
in HBM, we undertook a descriptive transcriptomic study of two 
HBM cases for which we had access to primary osteoblast cultures. 
When expression levels were compared between the HBM cases 
and five controls, it was interesting to find four genes that 
presented mRNA levels displaying a negative trend with BMD 
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Figure 3. Analysis of mRNA levels of several candidate genes in relation to BMD levels. (A-D) Trend of correlation betw/een Z-score values 
and gene expression levels of (A) JWISTI, (B) IL6R, (C) DLX3 and (D) PPARG. (E) One of the HBM samples presented an expression level of S0X6 5-fold 
decreased in relation to the mean of five control individuals. (G) The other HBM sample presented an expression level of RUNX2 that was 6-fold 
increased. 

doi:10.1371/journal.pone.0094607.g003 
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{TWISTl, IL6R, DLX3 and PPARG). Two other genes had one of 
the two HBM samples with oudier mRNA levels: RUNX2 was 6- 
fold elevated in one HBM sample, while S0X6 was 5-fold reduced 
in the other HBM sample. In spite of the evident low sample size, 
it is tempting to speculate that these genes (among others) may 
play a role in HBM, acting in an additive way. 

The limitations of this study include the small sample size of 
HBM cases, which precluded reaching significance. Replication in 
other cohorts will be necessary to confirm some of our results. This 
is particularly important in the expression analysis, where the 
number of analysed individuals was modest because of the great 
difficulty in finding primary human osteoblasts. It may be argued 
that the fact that the donors suffered from arthritis is a 
confounding factor. However, all donors were affected with this 
condition regardless of their BMD and it might be assumed that 
the effect of artritis would be similar in all of the samples. To 
minimize the impact of arthritis, caution was taken to obtain the 
osteoblasts from locations far away from the lesion. However, a 
systemic effect of arthritis cannot be totally ruled out. 

In conclusion, to our knowledge, this genomic and transcrip- 
tomic analysis of HBM is the first report of its kind. By combining 
both strategies, it was possible to gain a deeper insight into the 
genetic makeup of HBM. It includes suggestive evidence of genetic 
heterogeneity based on the observation of additive effects of 
several genes, on one hand, and monogenic cases not caused by 
LRP5, on the other. DKKl, possibly responsible for one of these 
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